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ABSTRACT 

A. new two-fluid mathematical model for fully three-dimensional steady solidification under 
the influence of an arbitrary acceleration vector and with or without an arbitrary externally applied 
steady magnetic field have been formulated and integrated numerically. The model includes Joule 
heating and allows for separate temperature-dependent physical properties within the melt and the 
solid. Latent heat of phase change during melting/solidification was incorporated using an enthalpy 
method. Mushy region was automatically captured by varying viscosity orders of magnitude 
between liquidus and solidus temperature. Computational results were obtained for silicon melt 
solidification in a parallelepiped container cooled from above and from a side. The results confirm 
that the magnetic field has a profound influence on the solidifying melt flow field thus changing 
convective heat transfer through the boundaries and the amount and shape of the solid accrued. This 
suggests that development of a quick-response algorithm for active control of three-dimensional 
solidification is feasible since it would require low strength magnetic fields. 


INTRODUCTION 

Possible means of devising controlling mechanisms for solidification and bulk single 
crystal growth from a melt is of practical importance to many industrial processes [Kosovic, 
Dulikravich and Lee 1991; Dulikravich, Kosovic and Lee 1992]. We will demonstrate the 
quantified effects of a fully three-dimensional solidification control using magnetic fields. If there 
are no electric charges in the melt and no external electric field is applied, a magnetohydrodynamic 
(MHD) model [Stuetzer 1962] can be applied. MHD flows without solidification have already 
been numerically analyzed in three-dimensions [Ozoe and Okada 1989; Lee and Dulikravich 1991]. 
It was only recently that various researchers have demonstrated numerically the effects of the 
magnetic field in two-dimensional [Salcudean and Sabhapathy 1990; Dulikravich, Kosovic and Lee 
1991; Kosovic and Dulikravich 1991] and in three-dimensional [Dulikravich, Ahuja and Lee 
1993a, 1993b; Dulikravich and Ahuja 1993] solidification. 

Our objective is to use a single system of governing equations in the entire domain which 
could consist of the melt, mixture of the melt and the solid (mushy region), or the solid alone. In 
other words, we want to use a single system of partial differential equations that treats both the 
melt and the solid as liquids, while allowing each of the two liquids to have its own set of 
temperature-dependent physical properties (density, heat capacity, thermal conductivity, electric 
permitiyity, thermal expansion, magnetic permeability, etc.). A crucial difference between the two 
liquids is that the liquid which models the solid phase (T < T^^us) must be assigned extremely 

high viscosity. Hence, velocities interior to the regions occupied by the liquid that models the solid 
phase will be practically zero. Consequently, we will refer to this artificially extremely viscous 
liquid as "solid", while the pure melt (T > T}jq U j c j us ) will be called "liquid". 
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Mass fraction of the liquid at any point in the domain determines locally to what extent 
should physical properties of the liquid or the solid be taken into account Since latent heat 
released or absorbed per unit mass of the mushy region is proportional to the local volumetric 
fraction, f, occupied by the liquid in a particular computational cell, this ratio is often modeled 
[Voller and Swaminathan 1991] as 


f ^liquid ^ Q ~ ^solidus (1) 

^liquid + V solid 01iquidus " ® solidus 

where the exponent "n" is typically 0.2 < n < 5. Here, the non-dimensional temperature is defined 
as 0 = (T - To)/AT 0 where T is the temperature measured in degrees. Typically, either 

AT 0 = Tiiquidus - T solidus and T 0 = T solidus or AT o = T hot * T cold ^ T 0 = T cold- T ^ e * ateril 
heat , L, is released in the mushy region (where Tiiquidus >T> T solidus) in proportion to the 

fraction f of the liquid in the mixture. Liquid density, pi is assumed to vary linearly as a function 
of non-dimensional temperature 


Pi' = 


1 + 


9(Pl/P0l) 

30 


q (0 - 0 q) = 1 * a ol (0 * e 0 ) 


(2) 


with a similar expression for the solid phase. We will use subscripts "1” and "s" to designate liquid 
and solid phase, respectively. If subscript "o" designates reference values, then the non- 
dimensionalization can be performed as follows 
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with similar expressions for the solid phase. Here, v, g, H and x are vectors of velocity, gravity 

acceleration, magnetic field and spatial position vector, respectively. Similarly, |i, k, c and a are 
coefficients of viscosity, heat conductivity, heat capacity and thermal expansion, respectively. 
Hydrodynamic pressure, coefficients of electric conductivity and magnetic permeability are 

designated with p, a and y, respectively. In this work we assumed that a and y do not vary with 

temperature (Oj = aoi/a 0 , c s = aos/c 0 , yj = yoi/y 0 and y s = yo</Yo)- Since the reference values 
designated with the subscript "o" are arbitrary, the non-dimensional numbers can be defined as 
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P o “o Igol AT„ l\ 
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H t = Yo Hollo — 


Then, adopting an extended Boussinesq approximation [Gray and Giorgini 1976] to MHD flows, 
the non-dimensional Navier-Stokes equations for phase-changing mixtures of two liquids 
[Dulikravich and Ahuja 1993] become: 

Mass conservation for two-phase MHD flows 

V** v* = 0 (8) 

Linear momentum conservation for two-phase MHD flows with thermal buoyancy and 
magnetic force 

P f£ + f pjV-(y* v* + p* I) + (1-0 p* V*-(y* v* + p * I) 

* 2 

= f [V*-4j (V*v* + (V*y*) T » + p* <?** H *)* H * + °T rr g * )] 

Ke ' 1 P m Re z 1 Re z 

* 2 

+(l-f)[V*.(^| (V*v*+ (V*v*) T )) +p* (y* -^-2 (V*x H*)x H*+ a* g*)] (9) 

PjjjRe Re 

Energy conservation for incompressible two-phase MHD flows including Joule heating 


(f P! (c el 0) i0 + (1-0 P s (c«e) e ) pr + f P 1 * V-tc^e y‘) + (1-0 fl. V*-(c es 9 y*) 


^ * * 


= f ( R^ V ** (k l V * 0) + ^ Ht 2 E< 3 (V*x H*)-( V*x H*» 

CT| Pm Re 

* 2 

+ (1-0 fejp! V*-(k* V*0) + -A- H ‘ 2 Ec , (V*x H*)-(V*x H*)) 

o s P m Re 

Magnetic field transport equations for two-phase MHD flows 

9H* „* * * f/(a*y*) + (1 - f)/(o*Ys) * 2 * 

yV- v*x (v*x H*) = — — — — V* 2 H* 

* p m R e P 


Here, p — f pj ^ (1 * f) p^ , where f — 1 for 0 ^ ^liquidus f — 0 for 0 < ^solidus* Non- 
dimensional hydrostatic, hydrodynamic, and magnetic pressures were combined to give 
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( 12 ) 
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where <J) is the non-dimensional gravity potential defined as g = V <j) . We used an enthalpy 
method to formulate the equivalent specific heat coefficients in the liquid and the solid phases are 
* * 1 9f . * * 1 9f 


c el = C 1 


andc.. 
^te 90 es 


= c s - — , respectively. This expression allows latent heat to be 


released in the mushy region according to the empirical law given in equation (1). 

At the solid walls, the velocity components were set to zero and the wall pressure was ----- 
determined from the normal momentum equation. Wall temperatures were either specified or 
obtained from the specified wall heat fluxes and the points on the first grid layer off the walls. 
Magnetic field was either specified as uniform on a particular wall, or its normal derivative to the 
wall was set to zero. The governing system of eight partial differential equations (8-11) was 
transformed into a non-orthogonal curvilinear coordinate system compatible with a typical three- 
dimensional boundary-conforming structured computational grid. In such a way, the resulting = 
finite difference algorithm for an iterative integration of the system can be applied to relatively 
arbitrary three-dimensional configurations. Since the system is singular (its time-dependent term in 
the mass conservation equation is zero), it can be integrated simultaneously only after introducing 
an artificially time-dependent term [Chorin 1967] in the mass conservation. Consequently, such an 
"artificial compressibility" iterative process does not follow physical time, but rather an artificial 
time coordinate. As a result, intermediate solutions are not time-accurate pictures of the flow field, 
but the final converged steady solution is accurate since the artificial compressibility term variation 
with iterations then becomes zero. 

We are using an explicit four-step time-integration and central differencing in space. Since 
the magnetic field transport equations (11) are strongly parabolic (for the given velocity field), their 
allowable integration time step is much smaller than in the case of the Navier-Stokes equations (7- 
10). Consequently, we coded the three magnetic transport equations separately from the three- 
dimensional Navier-Stokes equations [Lee and Dulikravich 1991] so that we can use different time 
steps for the two systems. Communication between the two systems based on periodically 
updating source terms (thermal buoyancy and magnetic effect s) in the Navier-Stokes system. 


NUMERICAL RESULTS 

Based on this analytical model for a two-fluid MHD solidification, a fully three-dimensional 
MHD flow analysis computer program was developed. Numerical results from this code were 
compared with known three-dimensional MHD analytical solutions in the case of no heat transfer 
[Lee and Dulikravich 1991] and in the case of heat transfer but without solidification [Ozoe and 
Okada 1991]. In both cases the code proved to be highly accurate [Dulikravich, Ahuja and Lee 
1993a]. This code was then augmented to incorporate temperature-dependent physical properties of 
the melt and the solid phase and the effects of latent heat release with an adequate account of the 
mushy region. 

We decided to study the three-dimensional MHD effects on solidification by numerically 
analyzing an MHD solidification in a parallelepipedal closed container initially filled with molten 

silicon. The container was discretized with 40 x 20 x 20 grid cells that were clustered 
symmetrically towards all the walls (Figure 1). Gravity was assumed to act vertically downward in 
the positive z-direction. If not indicated otherwise, the solid walls were thermally insulated. The 
values of the reference parameters were: lv 0 l = 0.02342 m s' 1 , 1 Q = 0.02 m, lg 0 l = 9.81 m s' 2 . 
Physical properties for silicon were compiled from a number of references (Table 1) which lead to 
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the following values of the non-dimensional numbers: Gr = 2.89 x 10^, Re = Gr^ = 1702, Ec = 
2.589 x 10* 8 , Pr = 0.01161, Pr = 4.255 x 10' 6 and Ht = 837.3 Bo, where Bo is measured in 
Teslas. The exponent used in the model for latent heat release (1) was n = 5. Test cases were run 
with Ht = 0, 20, 30, 60 corresponding to the magnetic fields having strengths of 0 Tesla, 0.02389 
Tesla, 0.03583 Tesla and 0.07166 Tesla, respectively. 

Solidification From the Top 

In this case the container had dimensions 0.04 m x 0.02 m x 0.02 m. Top wall was 
uniformly cooled below freezing temperature (0 = - 0.5) and the bottom wall was uniformly heated 
(0 = 0.5). A solidification case was first run without the magnetic field (Ht = 0). The computed 

contours of equal vertical velocity magnitude evaluated in the z* = 0.3 horizontal plane (Fig. 2a) 
indicate strong upward melt motion mainly along the short vertical walls and a centrally located 

downward jet. This is more evident in Figures 2b (evaluated at y* = 0.5 vertical mid-plane) and 2c 

(evaluated at x* = 0.5 vertical mid-plane) where it is clear that almost one-third (4926 solidified 
cells out of 16000 computational cells in the container) solidified starting from the top wall [Kerr et 
al. 1990]. Evidently, heat transfer is carried out by both conduction and convection. The computed 

isotherms in the vertical y* = 0.5 mid-plane (Fig. 2d) and in the vertical x* = 0.5 mid-plane (Fig. 
2e) indicate that the solid/melt interface is somewhat pulled down in the central part of die container 
due to the strong centrally located downward melt jet After reaching the bottom of the container, 
the jet spreads out and starts mowing upwards along the side walls thus forming a deformed 
vertical thoroidal melt motion. 

A uniform steady magnetic field of Ht = 20 was then assigned in the vertically downward 
direction (same as the gravity direction). The computed contours of equal vertical velocity 

magnitude evaluated in the z* = 0.25 horizontal plane (Fig. 3a) now indicate that the peak upward 
melt motion is mainly along the long vertical walls, while the centrally located downward jet 
became narrower and developed a non-parabolic profile. It is especially important to notice that the 
velocity profiles close to the walls and the solid/melt interface became steeper (Fig. 3b and 3c) 
which is typical of MHD flows. This in turn caused enhanced heat convection in the mushy region 
resulting in less accrued solid (4773 solidified computational cells). Because of the stronger 
downward centrally located jet the computed isotherms in the vertical mid-planes indicate slight 
sagging of the solid/melt interface (Fig. 3d and 3e). 

Finally, a uniform steady magnetic field of Ht = 60 was assigned in the vertically 
downward direction. The computed contours of equal vertical velocity magnitude evaluated in the 

z* = 0.25 horizontal plane (Fig. 4a) indicate a dramatic change as compared to the two previous 
cases. Specifically, the peak upward melt motion is contained in a narrow region along the short 
vertical walls, while the centrally located downward jet became wider and developed a strongly 
double-parabolic profile (Fig. 4a). The magnitudes of velocity vector components in the entire flow 
field were substantially reduced (Fig. 4b and 4c) to about 30% of those in the case with Ht = 20. 
Because of the weaker and less concentrated downward centrally located jet the computed isotherms 
in the vertical mid-planes indicate that the solid/melt interface is more planar (Fig. 4d and 4e). In 
this case there were 4676 solidified computational cells. 

Solidification From a Side 

The solidification was then tested for the case where one vertical wall (x* = 0) was kept 
uniformly hot (0 = 0.5) and the opposite vertical wall (x* = 1) was at a uniformly below freezing 

temperature (0 = - 0.5). Three runs were performed; one without the magnetic field (Ht = 0) and 
the other two with a uniform external magnetic field (Ht = 30) applied in the vertical z-direction and 
horizontal x-direction, respectively. The container size in this case was 0.01 m x 0.01 m x 0.02 m. 
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In the case of no magnetic field (Ht = 0), the computed contours of equal z-momentum 
evaluated in the z* = 0.5 horizontal mid-plane indicate (Fig. 5a) a single large vortex that rises the 
melt at the hot wall and descends it at the solid/melt interface (Fig. 5b). The computed isotherms in 

the vertical y* = 0.5 (Fig. 5c) and x* = 0.5 (Fig. 5d) mid-planes show that the solid/melt interface 
is highly curved and three-dimensional. In this case there were 3990 solidified computational cells. 
When the uniform magnetic field of Ht = 30 was applied downward, the computed contours 

of equal z-momentum evaluated in the z* = 0.5 horizontal mid-plane indicate (Fig. 6a) that the 
intensity of the single large vortex decreased (Fig. 6b). Also, a recirculation mini-vortex that 
existed at the lower comer of the hot vertical wall in the case without a magnetic field (Fig. 5b) was 
now eliminated. Computed isotherms in the vertical mid-planes (Fig. 6c and 6d) indicate that the 
solid/melt surface in this case became essentially two-dimensional. The number of solidified 
computational cells in this case increased significantly to 4283. 

A uniform magnetic field of Ht = 30 was then applied horizontally in the hot-to-cold x- 

direction. The computed contours of equal z-momentum evaluated in t he z * = 0.5 horizontal mid- 
plane indicate (Fig. 7a) that the intensity of the single large vortex decreased even further (Fig. 7b). 
Computed isotherms in the vertical mid-planes (Fig. 7c and 7d) indicate that the solid/melt interface 
became essentially two-dimensional. The number of solidified computational cells in this case 
increased significantly to 43 1 1 . 
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pi [kg m'3] 

2550 

p s [kg m' 3 ] 

2330 

q [J kg'l K'l] 

1059 

c s [Jkg-lK'l] 

1038 

ki [W m'l K'l] 

64 

k s [W m'l K'l] 

22 

Ti[K] 

1685 

T m [K] 

1683 

T S [K] 

1681 

jj.[kgm'l s'l] 

7.018 x 10' 4 

oci [K-l) 

1.41 x 10' 4 

Cts [K- 1 ] 

1.41 x 10- 4 

L [J kg" 1 ] 

1803000 


ci [W'l m*l] 
c s [W*l m"l] 
Y [Tm A'l] 


12.3 x 10^ 

4.3 x 10 4 
47t x 10-6 


Table 1. Physical properties used for silicon. 
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Figure 1. Computational grid and coordinate system 





















Figure 2. Solidification from the top: Ht = 0. 

a) constant vertical velocity component contours at z* - 0.3 horizontal plane; b) velocity vector field 
in vertical y* = 0.5 mid-plane; c) velocity vector field in vertical x* = 0.5 mid-plane; d) isotherms in 
vertical y* = 0.5 mid-plane; e) isotherms in vertical x* = 0.5 mid-plane. 




Figure 3. Solidification from the top: Ht = 20. 

a) constant vertical velocity component contours at z = 0.3 horizontal plane; b) velocity vector field 
in vertical y* - 0.5 mid-plane; c) velocity vector field in vertical x* = 0.5 mid-plane; d) isotherms in 
vertical y* = 0.5 mid-plane; e) isotherms in vertical x* = 0.5 mid-plane. 






Figure 4. Solidification from the top: Ht = 60. 

a) constant vertical velocity component contours at z* = 0.3 horizontal plane; b) velocity vector field 
in vertical y = 0.5 mid-plane; c) velocity vector field in vertical x* = 0.5 mid-plane; d) isotherms in 
vertical y* = 0.5 mid-plane; e) isotherms in vertical x* = 0.5 mid-plane. 
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Figure 5. Solidification from the side: Ht * 0. 

a) constant vertical velocity component contours at z* = 0.5 horizontal mid-plane; b) velocity vector 
field in vertical y* = 0.5 mid-plane; c) isotherms in vertical y* = 0.5 mid-plane; d) isotherms in 
vertical x = 0.5 mid-plane. 



Figure 6. Solidification from the side: Ht = 30 applied vertically. 

a) constant vertical velocity component contours at z* - 0.5 horizontal mid-plane; b) velocity vector 
field in vertical y* ■ 0.5 mid-plane; c) isotherms in vertical y* = 0.5 mid-plane; d) isotherms in 
vertical x = 0.5 mid-plane. 




Figure 7. Solidification from the side: Ht = 30 applied in hot-to-cold direction, 
a) constant vertical velocity component contours at z* — 0.5 horizontal mid-plane; b) velocity vector 
field in vertical y* = 0.5 mid-plane; c) isotherms in vertical y* = 0.5 mid-plane; d) isotherms in 
vertical x = 0.5 mid-plane. 
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